Cream of the Crop 1/Cream of the Crop 1.iso/PROGRAM/TPL60N14.ARJ/MAINVARS.PAS

Pascal/Delphi Source File | 1992-02-27 | 18KB | 581 lines

{$a+,n-,x-,s-,i-,r-,b-,v-}
unit mainvars;

interface
{ (C) Apr 19 1983 in BASIC version by:
      Professor W M Kahan,
      567 Evans Hall.
      Electrical Engineering & Computer Science Dept.
      University of California
      Berkeley, California 94720
      USA
 converted to Pascal by:
      B A Wichmann
      National Physical Laboratory
      Teddington Middx
      TW11 OLW
      UK
 further massaging by dmg =
      David M. Gay
      AT&T Bell Labs
      600 Mountain Avenue
      Murray Hill, NJ 07974
 and a couple of bug fixes from dgh = sun!dhough (29 May 1986)

 See the article by Richard Karpinski in the February 1985 issue
 of BYTE Magazine.

 You may copy this program freely if you acknowledge its source.
 Comments on the Pascal version to NPL or dmg, please.  }

   const
      {integer constants}
      NoTrials = 20;
      {Number of tests for commutativity. }

   type
      Guard = (Yes, No);
      Rounding = (Chopped, Rounded, Other);
      Message = packed array [1..40] of char;
      WhichOp = packed array [1..14] of char;
      Class = (Flaw, Defect, SeriousDefect, Failure);

   var
      {input: text;}
      {Small floating point constants.}
      Zero, { 0.0; }
      Half, { 0.5; }
      One, { 1.0; }
      Two, { 2.0; }
      Three, { 3.0; }
      Four, { 4.0; }
      Five, { 5.0; }
      Eight, { 8.0; }
      Nine, { 9.0; }
      TwentySeven, { 27.0; }
      ThirtyTwo, { 32.0; }
      TwoForty, { 240.0; }
      MinusOne, { -1.0; }
      OneAndHalf: { 1.5; } real;
      MyZero: integer;
      NoTimes, Index: integer;
      ch: char;
      AInverse, A1: real;
      Radix, BInverse, RadixD2, BMinusU2: real;
      C, CInverse: real;
      D, FourD: real;
      E0, E1, Exp2, MinSqrtError: real;
      SqrtError, MaxSqrtError, E9: real;
      Third: real;
      F6, F9: real;
      H, HInverse: real;
      I: integer;
      StickyBit, J: real;
      M, N, N1: real;
      Precision: real;
      Q, Q9: real;
      R, R9: real;
      T, Underflow, S: real;
      OneUlp, UnderflowThreshold, U1, U2: real;
      V, V0, V9: real;
      W: real;
      X, X1, X2, X8, RandomNumber1: real;
      Y, Y1, Y2, RandomNumber2: real;
      Z, PseudoZero, Z1, Z2, Z9: real;
      NoErrors: array [Class] of integer;
      Milestone: integer;
      PageNo: integer;
      GMult, GDiv, GAddSub: Guard;
      RMult, RDiv, RAddSub, RSqrt: Rounding;
      Continue, Break, Done, NotMonot, Monot, AnomolousArithmetic, IEEE,
            SquareRootWrong, UnderflowNotGradual: Boolean;
      { Computed constants. }
      {U1  gap below 1.0, i.e, 1.0-U1 is next number below 1.0 }
      {U2  gap above 1.0, i.e, 1.0+U2 is next number above 1.0 }

   procedure Page;
   function Int (X: real): real;
   function Sign (X: real): real;
   procedure Pause;
   procedure Instructions;
   procedure Heading;
   procedure Characteristics;
   procedure History;
   procedure notify(T: WhichOp);
   procedure TestCondition (K: Class; Valid: Boolean; T: Message);
   function Random: real;
   procedure SqrtXMinX (ErrorKind: Class);
   procedure NewD;
   procedure SubRout3750;
   function Power (X, Y: real): real;
   procedure DoesYequalX;
   procedure SubRout3980;
   procedure PrintIfNPositive;
   procedure TestPartialUnderflow;

implementation

   procedure Page;
     begin
       (* write(#$C) *) {FF in TURBO Pascal} writeln; writeln;
       end;

   function Int (X: real): real;
   {  simulates BASIC INT-function, which is defined as:
      INT(X) is the greatest integer value less than or
      equal to X. }


      function LargeTrunc (X: real): real;

         var
            start, acc, y, p: real;
            trunced: integer; (* dgh *)

         begin (* LargeTrunc *)
         if abs (X) < maxint then begin
            trunced := trunc(X);
            LargeTrunc := trunced;
            end
         else
            begin
            start := abs (X);
            acc := 0.0;
            repeat
               y := start;
               p := 1.0;
               while y > maxint - 1.0 do
                  begin
                  y := y / Radix;
                  p := p * Radix;
                  end;
               trunced := trunc(y);
               acc := acc + trunced * p;
               start := start - trunced * p;
            until start < 1.0;
            if X < 0.0 then
               LargeTrunc := - acc
            else
               LargeTrunc := acc
            end;
         end (* LargeTrunc *);


      begin (* Int *)
      if X > 0.0 then
         Int := LargeTrunc (X)
      else if LargeTrunc (X - 0.5) = X then
         Int := X
      else
         Int := LargeTrunc (X) - 1;
      end (* Int *);


   function Sign (X: real): real;

      begin (* Sign *)
      if X < 0.0 then
         Sign := - 1.0
      else
         Sign := + 1.0;
      end (* Sign *);


   procedure Pause;

      var
         ch: char;

      begin (* Pause *)
      writeln ('To continue, press any key and newline:');
      readln (input);
      while not eoln (input) do
         read (input, ch);
      Page;
      write ('Diagnosis resumes after milestone no ', Milestone);
      writeln ('               Page: ', PageNo);
      writeln;
      Milestone := Milestone + 1;
      PageNo := PageNo + 1;
      end (* Pause *);


   procedure Instructions;

      begin (* Instructions *)
      writeln ('Lest this program stop prematurely, ',
            'i.e. before displaying');
      writeln ('         "END OF TEST",');
      writeln ('try to persuade the computer NOT to',
            ' terminate execution whenever an');
      writeln ('error like Over/Underflow or Division by Zero occurs,',
            ' but rather');
      writeln ('to persevere with a surrogate value after, ',
            ' perhaps, displaying some');
      writeln ('warning.  If persuasion avails naught, don''t despair'
            , ' but run this');
      writeln ('program anyway to see how many milestones it passes,',
            ' and then');
      writeln ('amend it to make further progress.');
      writeln ('Answer questions with Y, y, N or n',
            ' (unless otherwise indicated).');
      writeln;
      end (* Instructions *);


   procedure Heading;

      begin (* Heading *)
      writeln ('Users are invited to help debug and augment',
            ' this program so it will');
      writeln ('cope with unanticipated and newly uncovered',
            ' arithmetic pathologies.');
      writeln ('Please send suggestions and interesting results to');
      writeln('        Richard Karpinski');
      writeln('        Computer Center U-76');
      writeln('        University of California');
      writeln('        San Francisco, CA 94143-0704, USA');
      writeln;
      writeln('In doing so, please include the following information:');
      writeln('        Version:  10 February 1989');
      writeln('        Computer:'); writeln;
      writeln('        Compiler:'); writeln;
      writeln('        Optimization level:'); writeln;
      writeln('        Other relevant compiler options:'); writeln;
      end (* Heading *);


   procedure Characteristics;

      begin (* Characteristics *)
      writeln (
            'Running this program should reveal these characteristics');
      writeln ('  Radix = 1, 2, 4, 8, 10, 16, 100, 256, or ...');
      writeln ('  Precision = number of significant digits carried.');
      writeln ('  U2 = Radix/Radix^Precision = One Ulp (OneUlpnit in the');
      writeln ('    Last Place) of 1.000xxx .');
      writeln ('  U1 = 1/Radix^Precision = One Ulp of numbers',
            ' a little less than 1.0 .');
      writeln ('  Adequacy of guard digits for Mult., Div., and Subt.');
      writeln ('  Whether arithmetic is chopped, correctly rounded, ',
            'or something else');
      writeln ('    for Mult., Div., Add/Subt. and Sqrt.');
      writeln ('  Whether a Sticky Bit is used correctly for rounding.');
      writeln ('  UnderflowThreshold = an Underflow Threshold.');
      writeln ('  E0 and PseudoZero tell whether underflow is abrupt,',
            ' gradual, or fuzzy.');
      writeln ('  V = an overflow threshold, roughly.');
      writeln ('  V0  tells, roughly, whether Infinity is represented.');
      writeln ('  Comparisions are checked for consistency with',
            ' subtraction');
      writeln ('    and for contamination with pseudo-zeros.');
      writeln ('  Sqrt is tested.  Y^X is not tested.');
      writeln ('  Extra-precise subexpressions are revealed',
            ' but NOT YET tested.');
      writeln ('  Decimal-Binary conversion is NOT YET tested',
            ' for accuracy.');
      end (* Characteristics *);


   procedure History;

      begin (* History *)
      writeln ('The program attempts to discriminate among');
      writeln ('   FLAWs, like lack of a sticky bit,');
      writeln ('   SERIOUS DEFECTs, like lack of a guard digit, and');
      writeln ('   FAILUREs, like 2+2 = 5 .');
      writeln ('Failures may confound subsequent diagnoses.');
      writeln;
      writeln ('The diagnostic capabilities of this program go beyond',
            ' an earlier');
      writeln ('program called "MACHAR", which can be found at the',
            ' end of the');
      writeln ('book  "Software Manual for the Elementary Functions',
            '" (1980) by');
      writeln ('W. J. Cody and W. Waite.  Although both programs',
            ' try to discover');
      writeln ('the Radix, Precision and range (over/underflow',
            ' thresholds)');
      writeln ('of the arithmetic, this program tries to cope',
            ' with a wider variety');
      writeln ('of pathologies, and to say how well the',
            ' arithmetic is implemented.');
      writeln ('BASIC version of this program (C) 1983 by Professor',
            ' W. M. Kahan;');
      writeln ('see source comments for more history.');
      writeln;
      writeln ('The program is based upon a conventional',
            ' radix representation for');
      writeln ('floating-point numbers, but also allows',
            ' logarithmic encoding');
      writeln ('as used by certain early WANG machines.');
      writeln;
      end (* History *);


   procedure notify(T: WhichOp);
      begin
      writeln('The ', T, ' test is inconsistent;');
      writeln('   PLEASE NOTIFY KARPINSKI !');
      end;

   procedure TestCondition (K: Class; Valid: Boolean; T: Message);

      begin (* TestCondition *)
      if not Valid then
         begin
         NoErrors [K] := NoErrors [K] + 1;
         case K of
            Flaw:
               write ('FLAW');
            Defect:
               write ('DEFECT');
            SeriousDefect:
               write ('SERIOUS DEFECT');
            Failure:
               write ('FAILURE');
            end;
         writeln (':  ', T);
         end;
      end (* TestCondition *);


   function Random: real;

      var
         X, Y: real;

      begin (* Random *)
      X := RandomNumber1 + R9;
      Y := X * X;
      Y := Y * Y;
      X := X * Y;
      Y := X - Int (X);
      RandomNumber1 := Y + X * 0.000005;
      Random := RandomNumber1;
      end (* Random *);


   procedure SqrtXMinX (ErrorKind: Class);

      begin (* SqrtXMinX *)
      SqrtError := ((sqrt (X * X) - X * BInverse) - (X - X * BInverse))
            / OneUlp;
      if SqrtError <> 0 then
         begin
         if SqrtError < MinSqrtError then
            MinSqrtError := SqrtError;
         if SqrtError > MaxSqrtError then
            MaxSqrtError := SqrtError;
         J := J + 1;
         writeln;
         if ErrorKind = SeriousDefect then
            write ('SERIOUS ');
         writeln ('DEFECT:  sqrt( ', X * X, ' - ', X, ') = ');
         writeln (OneUlp * SqrtError, ' instead of correct value 0 .');
         end;
      end (* SqrtXMinX *);


   procedure NewD;

      begin (* NewD *)
      X := Z1 * Q;
      X := Int (Half - X / Radix) * Radix + X;
      Q := (Q - X * Z) / Radix + X * X * (D / Radix);
      Z := Z - Two * X * D;
      if Z <= Zero then
         begin
         Z := - Z;
         Z1 := - Z1;
         end;
      D := Radix * D;
      end (* NewD *);


   procedure SubRout3750;

      begin (* SubRout3750 *)
      if not ((X - Radix < Z2 - Radix) or (X - Z2 > W - Z2)) then
         begin
         I := I + 1;
         X2 := sqrt (X * D);
         Y2 := (X2 - Z2) - (Y - Z2);
         X2 := X8 / (Y - Half);
         X2 := X2 - Half * X2 * X2;
         SqrtError := (Y2 + Half) + (Half - X2);
         if SqrtError < MinSqrtError then
            MinSqrtError := SqrtError;
         SqrtError := Y2 - X2;
         if SqrtError > MaxSqrtError then
            MaxSqrtError := SqrtError;
         end;
      end (* SubRout3750 *);

   { Sun version of Power (from dgh)...
   function pow(x, y : real ): real; external c;
   function Power(X, Y: real): real;
      var xx, yy: real;
      begin
         xx := X;
         yy := Y;
         Power := pow(xx,yy);
         end;
   }

   function Power (X, Y: real): real;

  var z:real;
    iy:integer;
    reciproc:boolean;
      (* Crude power function:  see Cody & Waite for a better one. *)

      begin (* Power *)
      if Y = 0.0 then
         Power := 1.0
      else if (X = 0.0) and (Y > 0.0) then
         Power := 0.0
      else if (Y = Int (Y)) AND (Abs(Y) <= MAXINT) THEN BEGIN
         IY := TRUNC (Y);
         reciproc := iy < 0;
         iy := abs (iy);
         z := 1.0;
         while iy > 0 do begin
           if odd(iy) then begin
             z := z*x;

           end else begin
           end; {endif}
           iy := iy div 2;
           if iy >0 then begin
              x := sqr (x);
           end else begin
           end; {endif}

         end; {endwhile}
         if reciproc then begin
           power := 1 / z;

         end else begin
           power := z;
         end {endif}
         end
      else if (X < 0.0) and (Y = Int(Y)) then
         if odd(trunc(Y)) then Power := -exp (Y * ln (-X))
         else Power := exp (Y * ln (-X))
      else
         Power := exp (Y * ln (X));
      end (* Power *);


   procedure DoesYequalX;

      begin (* DoesYequalX *)
      if Y <> X then
         begin
         if N <= 0 then
            begin
            if (Z = Zero) and (Q <= Zero) then write('WARNING') (* dgh: Y --> Z *)
            else begin
               NoErrors [Defect] := NoErrors [Defect] + 1;
               write('DEFECT');
               end;
            writeln (':  computed (', Z, ') ^ (', Q, ') = ');
            writeln (Y, ', which compares unequal to correct ', X, ';'); (* dgh: V --> Y *)
            writeln (' they differ by ', Y - X);
            end;
         N := N + 1;
      { ... count discrepancies. }
         end;
      end (* DoesYequalX *);


   procedure SubRout3980;

      begin (* SubRout3980 *)
      repeat
         Q := I;
         Y := Power (Z, Q);
         DoesYequalX;
         I := I + 1;
         if I <= M then
            X := Z * X;
      until (X >= W) or (I > M);
      end (* SubRout3980 *);


   procedure PrintIfNPositive;

      begin (* PrintIfNPositive *)
      if N > 0 then
         writeln ('Similar discrepancies have occurred ', N, ' times.');
      end (* PrintIfNPositive *);


   procedure TestPartialUnderflow;

      begin (* TestPartialUnderflow *)
      N := 0;
      if Z <> 0 then
         begin
         writeln ('Since comparison denies Z = 0, evaluating');
         writeln ('(Z + Z) / Z should be safe.');
         write ('What the machine gets for (Z + Z) / Z is: ');
         Q9 := (Z + Z) / Z;
         writeln (Q9);
         if (abs (Q9 - Two) < Radix * U2) then
            begin
            write ('This is O.K., provided Over/Underflow');
            writeln (' has NOT just been signaled.');
            end
         else if (Q9 < One) or (Q9 > Two) then
            begin
            N := 1;
            NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
            writeln ('This is a VERY SERIOUS DEFECT!');
            end
         else
            begin
            N := 1;
            NoErrors [Defect] := NoErrors [Defect] + 1;
            writeln ('This is a DEFECT.');
            end;
         V9 := Z * One;
         RandomNumber1 := V9;
         V9 := One * Z;
         RandomNumber2 := V9;
         V9 := Z / One;
         if (Z = RandomNumber1) and (Z = RandomNumber2)
               and (Z = V9) then
            begin
            if N > 0 then
               Pause
            end
         else
            begin
            N := 1;
            NoErrors [Defect] := NoErrors [Defect] + 1;
            writeln ('DEFECT: What prints as Z = ', Z, 'compares');
            write ('different from        ');
            if not (Z = RandomNumber1) then
               writeln ('Z * 1 = ', RandomNumber1);
            if not ((Z = RandomNumber2)
                  or (RandomNumber2 = RandomNumber1)) then
               writeln ('1 * Z = ', RandomNumber2);
            if not (Z = V9) then
               writeln ('Z / 1 = ', V9);
            if RandomNumber2 <> RandomNumber1 then
               begin
               NoErrors [Defect] := NoErrors [Defect] + 1;
               writeln ('DEFECT   Multiplication does not commute;');
               writeln ('comparison alleges that 1 * Z = ',
                     RandomNumber2);
               writeln ('differs from Z * 1 = ', RandomNumber1);
               end;
            if N > 0 then
               Pause;
            end;
         end;
      end (* TestPartialUnderflow *);
end.